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ABSTRACT 

In experimental microwave maps, point-sources can strongly affect the estimation of the power-spectrum and/or the 
test of Gaussianity of the Cosmic Microwave Background (CMB) component. As a consequence, their removal from the 
sky maps represents a critical step in the analysis of the CMB data. Before removing a source, however, it is necessary 
to detect it and source extraction consists of a delicate preliminary operation. In the literature, various techniques have 
been presented to detect point-sources in the sky maps. The most sophisticated ones exploit the multi-frequency nature 
of the observations that is typical of the CMB experiments. These techniques have "optimal" theoretical properties 
and, at least in principle, are capable of remarkable performances. Actually, they are rather difficult to use and this 
deteriorates the quality of the obtainable results. In this paper, we present a new technique, the weighted matched filter 
(WMF), that is quite simple to use and hence more robust in practical applications. Such technique shows particular 
efficiency in the detection of sources whose spectra have a slope different from zero. We apply this method to three 
Southern Hemisphere sky regions - each with an area of 400 deg 2 - of the seven years Wilkinson Microwave Anisotropy 
Probe (WMAP) maps and compare the resulting sources with those of the two seven-year WMAP point-sources 
catalogues. In these selected regions we find seven additional sources not previously listed in WMAP catalogues and 
discuss their most likely identification and spectral properties. 

Key words. Methods: data analysis - Methods: statistical - Cosmology: cosmic microwave background 



i 1. Introduction most promising ones. A good example is the multi-frequency 

' matched filter (MMF), a well known technique in the com 

— I The detection of point-sources embedded in a noise back- munity of the 11 digital signal processing" (e.g. see |Ka^ 



ground is a critical issue in the analysis of the experimen- UM> that has been recently proposed by | Herranz et al 
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tal Cosmic Microwave Background (CMB) maps. The es- d^ffil and| Lanz et al. | d2Qlfl). Although in principle such 
timation of the power-spectrum of the CMB component technique has "optimaF properties, its use is rather difficult 



H ' and the test of its possible non-Gaussian nature need "a limiting the actual performance m real experimental scenar- 

9 s . ! priori" detection and removal of these sources. In partic- 10S - In addition, most of the detection methods available m 

ular the former operation is rather delicate due to the literature have been developed in the context of full-sky 

usually present diffuse background of astrophysical nature observations. These experiments undoubtedly represent an 

and/or the inevitable instrumental noise. With the increas- important tool to better understand the physical properties 

ing high sensitivities of instruments to detect the CMB of CMB - However, they suffer the drawback of the Galactic 

signal, the astrophysical foregrounds have become the ma- contamination that, in spite of the optimism expressed in 

jor source of contamination with respect to the instru- man y Papers, it could not been yet completely removed. In 

mental noise. In CMB experiments, the foreground sig- this work we present a new technique, that we call weighted 

nals at high galactic latitudes come mainly from the emis- snatched filter (WMF), tailored to the detection of point- 

sion of extragalactic point-sources. Given its importance, sources m regions where the Galactic contamination can 

this sub ject has been extensively considered in literature be considered negligible. This technique takes also into ac- 

(see e.g. lHerranz and Sanz ll2008at ICaravalho et al. II2009L count the instrumental noise and can be applied to small 
and references therein). Among the various proposed tech- regions, 
niques, the multi-frequency approaches appear to be the 
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Observational radio data such as those provided by the 
Wilkinson Microwave Anisotropy Probe (WMAP) satellite 
should be a test to the robustness of the WMF to de- 
tect extragalactic point-sources. Since this technique can be 
used to small sky patches, we plan in the future to exploit 
the high sensitivity and angular resolution of the Atacama 
Large Millimeter/submillimeter Array (ALMA) and to test 
its applicability to this facility we plan to apply the WMF 
to simulated ALMA data (Ramos et al., in preparation). 

The paper is divided as follows: in Sec. [5] the mathemat- 
ics of the WMF is described. Associated numerical experi- 
ments are shown in Sec. [3l The application of the WMF to 
the WMAP observational data and the identification of the 
new discovered sources are presented in Sec.|H Conclusions 
and future development of this work are discussed in Sec. 
[5j To speed up the reading of the article most of the math- 
ematical technical details are deferred to the Appendix. 



2. Point-source detection at high Galactic latitude 

In the context of point-source detection, data can be 
thought as two-dimensional discrete maps {Afj}-^ x , each 
of them containing N p pixels, corresponding to M different 
observing frequencies (channels), with the form 



Xi = Si + AT, 



(1) 



Here, <Si corresponds to the contribution of the point- 
sources at the ith frequency, whereas Afi denotes the cor- 
responding noise component. At high Galactic latitudes, 
the CMB component is expected to be the dominant one. 
Hence, Afi may be modeled by 



Afi = B + £ t , 



(2) 



where B is the contribution of the CMB component which 
is the same at all frequencies (in terms of thermodynamic 
temperature units) and £i is the instrumental noise corre- 
sponding to the ith channel. 

The contribution of the point-sources is assumed to have 
the form 

Si = diG, (3) 

with di the amplitude of the source to the ith channel and 
where all the sources are assumed to have the same profile 
Q independently of the observing frequency. Although, in 
general, this will not be true, it is possible to meet this con- 
dition by convolving the images with an appropriate kernel. 
In the following, the components {£i} are considered re- 
alizations of a Gaussian, stationary, zero-mean, stochastic 
process. 

The main feature of model CO)-© is that the CMB con- 
tribution does not change with the frequency. Hence, it is 
possible to linearly combine the maps in a single map 

y in such a way that the CMB contribution is zeroed. In 
particular, 

Al 
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WkXi, 



(4) 



where the weights w are chosen in order to fulfill the criteria 

(5) 
(6) 



w T l = 0, 



w T a = 1, 



with w = [u>i, u>2, . . . , wm] T , o, — [ai, d2, . . • , o,m} t , and 
1 = (1, 1, ....1) T . Here, symbol " T " denotes the matrix 
transpose. The first constraint ([3]) implies that the con- 
tribution of 3 in y is completely removed, whereas the 
second one © provides a normalizing factor. The advan- 
tage of this procedure is to deal with a map contaminated 
only by the instrumental noise. It is easier to deal with this 
kind of noise than with B since the correlation length of Af 
is much shorter than that of B. This is particularly useful 
in situations where only small patches of sky are available 
and the auto-covariance function of the noise (a piece of 
information necessary to any detection method) has to be 
estimated from the data. Moreover, working with a sin- 
gle map allows to use detection techniques as the classical 
matched filter (MF) (see Appendix whose robustness 
is proved by many years of applic ations in many d i fferen t 
fields of science and engineering (I Vio et al. I[200ll200l . 
We call here weighted matched filter (WMF) the coupling 
of the weighted combination of the maps with the MF. 

In the case of M — 2, (i.e. two maps are available), the 
only possible solution is w T = [l/(ai — a 2 ), — l/(ai — 02)]. 
However, for M > 2 more degrees of freedom are available. 
This allows the selection of the weights in such a way that 
specific conditions are satisfied. In particular, one could 
wish that the peak signal-to- noise ratio of y, 



R(w\a) 



(w T a) 2 
w T Dw ' 



is maximized, i.e. Q 



w = argmaxi?(ti;|a) 



(8) 



Here, D is the M x M cross-covariance matrix of the noise 
processes whose (i, j)th entry (D)ij is given by 



(D) 



0) 



with aZ; the variance of Si and af, the covariance between 

LI L J 

Si and £j . Because of the constraint (JB]) , condition can 
be reformulated as 



w = argmin[ti; T £)K;]. 



(10) 



This approac h diffe rs from that proposed by 
IChen and Wright! ([2009) which consists in the mini- 
mization of the simpler quantity w T w (i.e. instrumental 
noise is not taken into account). Moreover, these authors 
seem to adopt a numerical approach for such operation 
(no details are provided respect to this) . Actually, a simple 
analytic solution of problem (ITUl) . with the constraints ([5|)- 
(O, can be obtained by means of the Lagrange multipliers 
method, i.e. 



(D - 1 ! 



where 



7] = 


1 T D H- 


c = 


a T D 


= 





(11) 

(12) 
(13) 
(14) 



We recall that the functions "arg min F(x)" and 
"arg max provide the values of x for which the function 

F(x) has the smallest and greatest value, respectively. 
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This r esult is similar to that obta ined bv lRemazeilles et al. I 
(|20irfl and|H urier et aTTI (|201dt ) however, in a completely 
different context, in which the problem of interest is the sep- 
aration of CMB and Sunyaev-Zel'dovich effect through the 
internal linear composition appr oach. The main dif f erence 
wit h the solution (JTTT) i s that in iRemazeilles et al."l (|2010t ) 
and lHurier et al. 1 (|2010h the empirical covariance matrix of 
the observed maps Xi is used, instead of the matrix D in 

CD- 

3. Numerical experiments 

To test the performances of the WMF we have carried 
out some numerical experiments. We consider a scenario 
where four different observing frequencies are available. 
We make the simplifying assumption that all the chan- 
nels have the same point-spread function (PSF) which is 
a two-dimensional circular symmetric Gaussian normalized 
to have a peak value equal to one and with a dispersion set 
to three pixels. Here the CMB component B is simulated 
on a regular two-dimensional grid containing (101 x 101) 
pixels with size 3.52' x 3.52'. The instrumental noise £i is 
assumed to be a Gaussian white-noise process with variance 
equal to one in units of the standard deviation of the CMB 
signal. This scenario mimics that expected for the "Low- 
Frequency Instru ment" mounted on the PLANCK satellite 
(|Vio et al. 112003ft . The amplitudes {ai} of the point-sources 
are assumed to follow a power-law 

a-i = (jf^j °i ( 15 ) 

where z/£, i = 1,2,3,4 are the observing frequencies (i = 
1 -> 30 GHz, i = 2 -> 44 GHz, i = 3 -> 70 GHz, i = 4 -> 
100 GHz), ai — 0.5 (in units of the standard deviation of 
the CMB signal) and a is the source spectral index. The 
value of ai has been chosen to reproduce an experimental 
situation characterized by a rather low signal-to-noise ratio. 

Figure [1] shows the "probability of detection" , Pd , 
against the " probability of false alarm" , Pfa (Le. the proba- 
bility of a false detection), for the WMF (see Appendix [A"| . 
Four different values of a are considered, i.e. a — 
3, 1, 0.5, 0.05 (negative values of a provide similar results). 
For comparison we show also the results obtained with 
the classical multi-frequency matched filter (MMF), and 
those obtained using the WMF with the weights w = 
[p,p,..., —(M - l)p] T , p = 1/y^l - M) + (1 - Mf. This 
last method, that we name uniformly weighted matched fil- 
ter (UWMF), corresponds to a situation where only one 
signal is used to eliminate the component B, whereas the 
others are given an identical weight. Quantity p is fixed in 
such a way that w T l = and w T w = 1. In this exper- 
iment MMF and UWMF are used as, respectively, upper 
and lower limit for the results obtainable by WMF. This 
is because, under the conditions we are working with, no 
detection method can outperform MMF (see Appendix |A"|). 
On the other side, in general, it is not expected that UWMF 
achieves good performances, since the weights are com- 
puted without considering the noise level in the map as 
well the characteristics of the source spectra. In this kind 
of diagram, a method is superior to another one when, for 
a fixed Pfa, the corresponding Pd is greater. More specif- 
ically, the relationship Pd against Pfa should always be 
well above a 45° straight line (the dashed line in the figure) 



since this corresponds to a detection performance identical 
to that of flipping a coin, ignoring all the data. 

From Fig. [TJ it is evident that, when a > 1, the WMF 
and MMF have very similar performances. When 0.5 < a < 
1, the behavior of these both filters is still reasonably simi- 
lar. In the case of a ps (i.e. for sources with flat spectrum) 
the performance of the WMF becomes close to that of the 
UWMF. This happens because for a close to zero, not only 
the CMB has the same contribution at the various frequen- 
cies but also the intensity of the sources is constant. In this 
case, a simple alternative is to average the maps and then 
apply the classic MF to the resulting map. We call this 
method average matched filter (AMF). A point to stress 
is that, similarly to the MMF, also the AMF requires the 
knowledge of the auto-covariance matrix of the CMB. This 
last method, however, has the advantage over the MMF be- 
cause only one map has to be handled. For comparison, in 
Fig. [1] the performance of the AMF is also shown. 

From these considerations, it appears that an effective 
and the simplest procedure to detect less common sources 
(those with spectra different from flat) consists of using the 
WMF, optimised for different values of a. The reliability of 
such procedure is supported by Fig. [5] that shows the rela- 
tive decrease of the probability of detection, (P£ — Pd)/ Pd, 
against Pfa when the WMF is applied to a source whose 
true spectral index a is erroneously assumed to be a* . Here, 
Pd and P D are the probability of detection when the WMF 
is applied assuming the true and the wrong spectral index, 
respectively. The set of values [3, 1, 0.5, 0.05] is used for both 
a and a*. From this figure, it is evident that a remarkable 
decrease of the probability of detection is to be expected 
only if a is quite different from a* . 

4. Application to WMAP data 

4.1. WMAP maps 

The Wilkinson Microwave Anisotropy Probe (WMAP) 
satellite (Bennett et al. 2003a) was designed to produce 
microwave full-sky maps of the CMB radiation. With the 
aim of separate the CMB signal from foreground compo- 
nents, the maps were obtained at five different frequency 
bands, respectively centered at 23GHz (K band), 33 GHz 
(Ka band), 41 GHz (Q band), 61 GHz (V band) and 94 
GHz (W band). Although the limited angular resolution 
of WMAP, with a full width at half maximum (FWHM) 
of about 13 arcminutes (W band), it is presently the only 
one that can offer a millimeter wavelengths all-sky survey, 
providing an unique tool for the study of radio sources. 

The seven-year full sky temperature and polarization 
maps per frequency band, namely the Stokes I, Q and U 
parameters, are available from the LAMBDA websitfl The 
maps of the five frequency bands have different resolutions, 
from roughly 0.21° (W band) to about 0.82° (K band). 
With the goal of testing the proposal WMF method to 
identify radio point-sources we use the Stokes I (tempera- 
ture) co-added maps (combination of the individual differ- 
encing assemblies of a single frequency band) to a common 
1° FWHM Gaussian beam, from which the CMB dipole 
has been removed (for more details see Jarosik et al. 2010). 
These maps were generated as a nested HEALPbJ^ sky pro- 

2 http://lambda.gsfc.nasa.gov/ 

3 http://healpix.jpl.nasa.gov/ 
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Fig. 1. Probability of detection Pd against probability of false alarm Pfa for the detection methods in the numerical 
experiment described in Sec. [3] Here, the sources are assumed to have intensity at the i-th frequency given by 
o-i = (Vi/Vi) Q ai, with a\ = 0.5. Four values for a are considered, a — 3, 1,0.5,0.05. The instrumental noise is assumed 
to be a Gaussian, zero-mean, white-noise process whose standard deviation is equal to one in units of the standard 
deviation of the CMB signal. The results are shown for the different methods: the weighted matched filter (WMF), the 
multi-frequency matched filter (MMF), the uniformly weighted matched filter (UWMF) and the average matched filter 
(AMF). The MMF is used as benchmark since it has the best theoretical detection performance. The UWMF shows the 
worst possible results obtainable with the WMF approach. The AMF shows what results are obtainable when the maps 
are simply averaged. A method is superior to another one when, for a fixed Pfa, the corresponding Pd is greater (for a 
given method, the relationship Pd against Pfa should always be well above a 45° straight line, the dashed line in the 
figure). With increasing of a, the behavior of MMF and WMF becomes similar. For an a close to zero, the MMF and 
AMF show a very similar performance. 



jection (Gorski et al., 2005) with a resolution of Nside=512 
(corresponding to the label WMAP resolution of Res 9). 

4.2. Selection of the sky regions 

The drawback to use the CMB experiments for cosmo- 
logical studies is the foreground contamination from the 
Galaxy and extragalactic sources. In particular, the extra- 
galactic point-sources contaminate the CMB maps at fre- 
quencies below 60 GHz and at high frequencies the statisti- 
cal properties are still barely known. At high Galactic lat- 



itudes (|b|>15°) and for frequencies between 30 and 150 
GHz, the CMB signal dominates the Galactic one (e.g. 
Bennett et al. 2003b, Tegmark et al. 2000). Therefore, 
a strategic way to obtain more accurate cosmological in- 
formation is to observe at high Galactic latitudes where 
the foreground contamination is expected to be lower. We 
selected three particular sky regions, each with an area 
of 20°x20° centered at galactic longitude (1) and lati- 
tude (b) of (l,b) = (258.18° ,-46.33°), (l,b) = (252.07°,-38.78°) 
and (l,b) = (272.48°,-54.63°) denoted here by first, second 
and third sky region, respectively. The first coordinates 
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Fig. 2. Relative decrease of the probability of detection, (P£> — Pd)/Pd , against Pfa when the WMF is applied to a source 
whose true spectral index a is erroneously assumed to be a* . Here, Pd and Pq are the probability of detection when the 
WMF is applied assuming the true and the wrong spectral index, respectively. The set of values [3, 1, 0.5, 0.05] is used for 
both a and a*. 



were chosen because t his region is of particular interest t o 
the EBEx experiment (jReichborn-Kiennerud et al. 1 12010). 
towards which we plan follow-up observations with the 
Atacama Large Millimetre/submillimetre Array (ALMA). 
The centres of the other two regions are two point-sources 
identified in the first one. We extracted these three sky re- 
gions from the smoothed full-sky maps per frequency band 
and project them in squared maps of (512 x 512) pixels us- 
ing the HEALPix software. In Fig. [3] the resulting maps in 
all the frequency bands for the first region are shown. We 
applied the WMF method to the three regions mentioned 
above and the linear composition map for the first region 
can be seen in Fig. [3] To fix the detection threshold, we con- 
sider an approach based on the empirical probability density 
function (EPDF) of the values of the pixels of the linear 
composition maps after the application of the matched fil- 
ter. This is the typically procedure used for the detection of 
point-sources in the CMB context. Many authors set such 
threshold to five times the standard deviation of the pixel 



values in the map (5a level). This method, however, suf- 
fers the impact o f the point-sources themselves (e.g. see 
iLeach et al.1l2008l ). For this reason, we have adopted a dif- 
ferent approach that is much less dependent on the amount 
of the spurious contributions. In particular, we claim a de- 
tection when the fluxes in the corresponding pixels have val- 
ues above a threshold given by the 98% percentile computed 
over all the pixels in the respective sky region. The choice 
of this approach is forced by the difficulties in the computa- 
tion of the detection threshold starting from an "a priori" 
Pfa since we could not used the level of the pixel noise 
provided by the WMAP team (defined as a — a /^N o \, s , 
where an and iV b s are values taken from the LAMBDA 
website □). The reason is that, before their linear composi- 
tion, the maps have been manipulated to obtain a common 
spatial resolution as well as to convert them from spherical 
to rectangular coordinates (our codes are developed for the 



4 http:/ /lambda. gsfc.nasa.gov/product/map/dr4/ 
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Fig. 3. Squared WMAP maps per frequency band for the first sky region, centered at (l,b)=(258.18°,-46.33°) with an 
area of 20°x20°. The map at the lower right corner corresponds to the linear composition map obtained with the WMF. 



case of small patches of sky). Moreover, the noise in the 
WMAP maps is not spatially uniform. The EPDF for the 
first region can be seen in Fig. @] It is evident that the bulk 
of the pixel values is confined in a rather restricted range, 
say [—10, 15] in internal units of our codes with a standard 
deviation of about a\§ = 2.98 units. Since the CMB com- 
ponent is not present in the linear composition map, it is 
reasonably to assume that these values are due only to the 
noise. The EPDF presents a long tail up to a value of about 
132 units that is due to the presence of the point-sources. 
The 98% empirical percentile of the EPDF roughly corre- 
sponds to the 5<7i5 level. Similar values are found for the 
other regions. 



4.3. Identification of the WMAP Point-Sources 

The WMAP seven-year Point Source Catalogs contains in- 
formation on the point-sources in the five frequency bands 
from 23 to 94 GHz, based on data from the first 7 years of 
the WMAP sky survey from 10 Aug 2001 to 9 Aug 2008, 
inclusive. The WMAP team has produced two point-source 
catalogues using different methods for the identification of 
the sources, namely the Five-band search technique and the 
three-band CMB-free technique 0- The former catalogue 
contains 471 point-sources and it is complete to 2 Jy for 
regions of the sky away from the Galactic plane. The latter 
catalogue was built using the three frequency bands from 
41 to 94 GHz. This last method identifies 417 point-sources 
in a linear combination map for which the weights were ob- 



http:/ /heasarc. gsfc.nasa.gov/W3Browse/all/wmapptsrc. html 
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Fig. 4. Empirical probability density function (EPDF) of the values of the pixels for the linear composition map of the 
first sky region after the application of the matched filter. The EPDF is different from zero in the range [—11, 132]. 
Here, only the range [—10, 30] is shown since outside such interval the EPDF is close to zero. The vertical red and green 
lines show the detection threshold based on the 98th percentile and the 5<7i5 level, respectively. Here, cris is the standard 
deviation of the pixels with values less or equal to 15. The cyan line provides the Gaussian probability density function 
with zero mean and standard deviation given by (T15. 



tained such that the CMB contribution was removed and 
point-sources with flat-spectrum were preserved. 

In the three sky regions selected for this work all the 
sources found with the WMF were cross-checked with 
those found in both WMAP seven-year catalogues. All the 
sources found by the WMAP team for the selected regions 
are detected by WMF. However, the WMF allows us to find 
more sources that are not listed in none of these catalogues. 
Those new sources can be seen in the Fig. [5] and Fig. [51 We 
have labeled them as Source 1 and 2 in the first region and 
3, 4, 5, 6 and 7 in the third region. All the sources detected 
by the WMF in the second region are listed already in the 
WMAP catalogues, except two {Source 1 and 2) that are 
already identified in the first region. 

As our detection method is sensitive to the source spec- 
tral index, we investigated what kind of spectral type ra- 



dio source populations were detected by WMAP. The Five- 
band catalogue provides an estimate of the spectral index 
(a) defined by a power law of the form, F oc v~ a , where 
F is the flux density and v is the frequency. We made use 
of this information to study the distribution of the spec- 
tral indices. The a lies in the range [-2.1,1.3] and three 
main spectral classes can be identified: about 81% of the 
population have a fiat-spectrum (-0.5<a<0.5), 16% show 
an inverted-spectrum (a< —0.5) and the remaining 3% a 
steep-spectrum (a>0.5). 

We test different values for a with the WMF and the 
new sources are detected for spectral indices in the range 
obtained from the Five-band catalogue, except in the case 
of assuming an a close to 0. As said in Sec. [3] the WMF 
is more efficient in finding sources with spectral index dif- 
ferent from zero. For strictly flat-spectrum sources (a ~ 0) 
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the method erases not only the CMB component but also 
the sources having equal intensity in all the WMAP bands. 
Therefore, this method is optimized for those sources hav- 
ing a flux dependent on the frequency. 

4.3.1. Sources fluxes 

To identify the new detected sources we need to have an 
estimation of the fluxes measured by WMAP. To this aim 
we used the original seven-years maps and integrate the 
flux density at the source location within the WMAP beam, 
assumed to have a Gaussian profile with FWHM as given 
in Hinshaw et al. (2009). Conversion factors from mK to 
Jy have been derived by comparing our derived fluxes with 
those given in the WMAP catalogues. 

4.3.2. Cross identification of the new discovered WMAP 
Point-Sources 

To identify possible counterparts we cross-correlated the 
new WMAP sources with catalogues found in the NED 
database, AT20G Catalogue (Murphy et al, 2010) and 
NEWPS sources (Massardi et al, 2009). We checked all 
the sources with a radio counterpart within a selected ra- 
dius of 12' which corresponds to 3 times the mean position 
uncertainty of the WMAP satellite (Gold et al, 2010). 

We summarize in Tab. [1] a list of possible counter- 
parts. The Fig. [7] and Fig. [5] show the reconstructed ra- 
dio Spectral Energy Distributions (SED) of the most likely 
counterparts with the WMAP data. The linear regression 
lines and respective spectral indices were obtained through 
the BCES(Y|X) ordinary least-squ ares method which takes 
into account measurement errors (Akritas and Bcrshadv. 
[1991 . 

Here we discuss briefly the most likely identifications of 
the new sources, which we associate to the brightest radio 
sources within our search radius. 

— Source # 1 The strongest radio sources within our 
search radius are MRC 0427-539B of which only one 
flux value at 408 MHz is available and IC 2082 that is 
a Galaxy pair (Gpair) with radio data from 408 MHz 
to 22 GHz. We plot the corresponding radio SED of 
both sources assuming that the WMAP fluxes belong 
to them. Both spectra are plotted in Fig. [7J In both 
cases we plot for reference also the best fit through the 
data which indicates a spectral behavior of an inverted- 
spectrum with a spectral index of about -0.60±0.07 and 
-1.00±0.09, respectively. 

— Source # 2 The most likely association is the radio 
source PKS 0437-454 which has data between 2.7 and 
150 GHz. Our reconstructed SED is reported in Fig. [7] 
Here the distribution of fluxes is complex and may show 
a SED with different components. The uncertainties re- 
lated to the identification and the source fluxes do not 
allow to make any further investigation with the present 
data. 

— Source # 3 The closest radio source with the strongest 
fluxes is PKS 0212-620 which is a candidate QSO with 
radio data from 843 MHz to 20 GHz. The fluxes are 
shown in Fig. [7] The source seems to be variable and 
differenc es in flux values at t he same frequency confirm 
this (see lSadler et al.. 1120061) . 



— Source # 4 The most likely association is PKS 0313- 
660 which has data from 843 MHz to 20 GHz and it is 
identified as a QSO. The SED built with WMAP data 
is reported in Fig. [7] and the fit gives the spectral index 
of ~ 0.15±0.01. 

— Source # 5 The most probable identification is the QSO 
PKS 0235-618 that has radio data in the range 408 MHz 
to 20 GHz. The SED built with WMAP data is reported 
in Fig. [7] resulting in a fit with a spectral index of ~ 
0.09±0.01. 

— Source # 6 There are three sources with bright radio 
fluxes, namely SUMSS J032356-602410, PKS 0322-605 
and PMN J0323-6026. We report the radio SEDs in 
Fig. [51 For the first and second possible counterpart, 
only one flux value is available. The PMN J0323-6026 
has data from 843 MHz to 2 GHz and is a variable 
QSO (see lSadler etaE1l2006l) . 

— Source # 7 The closest and most likely association is 
PKS 0226-559 with radio data from 843 MHz to 8.4 
GHz identified as a flat s pectrum (at freque ncies smaller 
than 8GHz) radio QSO (jHealev et al.ll2007[ ). The recon- 
structed SED with the WMAP data is shown in Fig. [H 
At high frequencies the spectrum appears to decrease 
with frequency and it is no longer flat. 

Most of all the other radio possible counterparts for which 
we do not give any SED (see Tab. [T]), have only one radio 
detection mostly at 843 GHz of the order of mJy or do not 
have published radio data. 

5. Conclusions 

In this paper we have presented a new technique, the 
weighted matched filter (WMF), to extract point-sources 
from astrophysical maps. This method is quite simple to 
use and more robust in practical applications and it is op- 
timal in extracting sources with a spectrum different from 
a flat one (i.e. with a spectral index close to zero). We have 
shown the reliability of this technique with some numerical 
simulations. 

We have then applied the method to three Southern 
Hemisphere sky regions - each with an area of 400 deg 2 
- of the seven-year WMAP temperature maps and com- 
pared the resulting sources with those of the two seven-year 
WMAP point-sources catalogues. 

We have found in these three regions seven additional 
sources not previously listed in WMAP catalogues and dis- 
cuss their most likely identification and spectral properties. 

We plan to investigate and explore further the applica- 
tion of the WMF technique and the identification of the 
new sources in future experiments, namely with Planck 
observational data and with simulations aimed at reproduc- 
ing the sky at the ALMA frequencies and spatial resolution. 
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Fig. 6. Same as Fig. [5] for the third sky region centered at (l,b)=(272.48°,-54.63°). 
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Table 1. New detected sources and possible counterparts. First column corresponds to the number of the new sources 
identified by arrows in Fig. [5] and Fig. |5] with the galactic coordinates (in degrees) in the second and third column. The 
fourth, fifth and sixth columns are the object name (ID) and galactic coordinates (in degrees) of the possible counterparts 
found in the NED database within a search radius of 12 arcminutes. The ID of the possible counterparts in the AT20G 
and NEWPS 5yr 3-sigma surveys is also present. The seventh column gives the distance (in arcminutes) between the 
WMAP coordinates and those of the possible counterparts and the last column gives the type of object. 
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SOURCE #1 (MRC 0427-539B + WMAP) 
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Fig. 7. Spectral Energy Distribution for the new sources using the derived WMAP data and NED data of the possible 
counterparts. The arrows represent upper limits corresponding to 3a, where a was obtained from the pixel noise. The 
WMAP data is plotted with data from MRC 0427-539B and IC 2082 for Source # 1, PKS 0437-454 for Source # 2, PKS 
0212-620 for Source # 3, PKS 0313-66019.0 for Source # 4 and PKS 0235-618 for Source # 5. For the Sources # 1, # 
4 and # 5 , it is also plotted the BCES(Y|X) ordinary least-squares regression lines. 
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Fig. 8. The same as Fig. [7] with WMAP data plotted with data from SUMSS J032356-602410, PKS 0322-605 and PMN 
J0323-6026 for Source # 6 and PKS 0226-559 for Source # 7. 
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Appendix A: The classic matched filter method 
to detect point-sources in CMB maps 



This appendix integrates the missing information 
of Sects. 2 and 3. We present here the classic matched 
filter method to detect point-sources in the case of 
multi-frequency observations. The goal is to facilitate 
the comparison of such approach with that proposed in 
this work. Arguments will be developed starting from 
the one-dimensional signal and single frequency case. 
The two-dimensional signal and multi-frequency case is 
developed in the second part. 

A.l. One-dimensional signal and single-frequency case 

In the case of the CMB observations, the following condi- 
tions are commonly assumed: 

1. The point-sources have a known spatial profile s — ag. 
The amplitude "a" is a scalar quantity different from 
source to source, whereas g is a function which, 
due to the instrument beam, is identical for all the 
sources. Function g is normalized in such a way that 
max{<7[0], g[l], . . . , g[N — 1]} = 1, where N is length of 
the function; 

2. The point-sources are embedded in a noise-background 
n, i.e. the observed signal x is given by x = s + n. In 
other words, noise is additive; 

3. Noise n is the realization of a stationary stochastic pro- 
cess with known covariance matrix 



C = E[nn T ]. 



Because of the Galactic contribution, especially at low 
Galactic latitudes, this hypothesis is not always satis- 
fied but we assume that it holds locally. This allows the 
computation of statistics such as the mean or the co- 
variance matrix that, otherwise, should not be possible. 
Without loss of generality, it is assumed that E[n] = 0. 

Under these conditions, the detection problem consists in 
deciding whether a; is a pure noise n (hypothesis Hq) or 
it contains also the contribution of a source s (hypothesis 
Hi). In this way, the source detection problem is equivalent 
to a decision problem where two hypotheses hold: 



Ho 
•Hi 



x 

X 



n 



(A.2) 



Under Ho the probability density function of x is given by 
p(x\Ho) whereas under Hi by p(x\Hi). At this point, it is 
necessary to fix the criterion to use for the detection, which 
depends on the particular case of interest. For example, one 
could decide that the non-detection or the misidentification 
of a bright source could be more important than the detec- 
tion of a fainter one, or viceversa. A very common and 
effective criterion is the Neyman-Pearson criterion which 
consists in the maximization of the probability of detection 
Pd under the constraint that the probability of false alarm 
Pfa (i.e., the probability of a false detection) does not ex- 
ceed ai^edvalue a. The Neyman-Pearson theorem (e.g., 
see [Ka^lH998h is a powerful tool that allows to design a 



decision process that pursues this aim: to maximize Pd for 
a given Pfa = ot, decide Hi if the likelihood ratio (LR) 



p(x\Hi) 
p(x\Hq) 

where the threshold 7 is found from 



Pfa = I p(x\Ho)dx = a. 

I {x:L(x)>~/} 



(A.3) 



(A.4) 



The test of the ratio (|A.3I) is called the likelihood ratio test 
(LRT). 

An important example of application of LRT is the case 
of a Gaussian noise n with correlation function C. Actually, 
in CMB experiments this condition is satisfied only for ob- 
servations at high Galactic latitudes where the CMB emis- 
sion and the instrumental noise are the dominant contribu- 
tions. At lower latitudes, it is often assumed to hold locally. 
For example, the contribution to x of components that in 
small sky patches show linear spatial trends are often ap- 
proximated with stationary Gaussian processes with a steep 
spectrum (e.g. 1// noises). This is usually assumed, even 
in cases of unrealistic Gaussianity condition, since it allows 
an analytical treatment of the problem of interest and the 
results can be used as a benchmark in the analysis of more 
complex scenarios. When n is Gaussian, 



p(x\Ho) = Aexp 
p(x\Hi) = Aexp 



l -x T C- x x 
2 



-{x - s) T C- 1 {x - s) 



(A.l) with 



A 



(27r)f det*(C)' 



(A.5) 
(A.6) 

(A.7) 



The LRT is given by 

l(x) = HL{x)] = x T C^ 1 s - is T C -1 s > i '. (A.8) 

Hence, it results that Hi has to be chosen when the statistic 
T(x) (called NP detector) is 



T(x) = x T C' 1 s > 7" 



with 7" such that 



Pfa = Q 



[s T C- x s] 



1/2 



(A.9) 



a, (A.10) 



i.e. 



7 



= Q- 1 (PfaWs t C- 1 s. 



(A.ll) 

Here, Q{x) is the complementary cumulative distribution 
function 



Q(x) 



1 I 1 2 ! , 

exp — t at, 
2tt V 2 1 



(A.12) 



and Q 1 the corresponding inverse function. 
Equation (|A.10[) is due to the fact that T(x) is a 
Gaussian random variable with variance s T C~ 1 s and 
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expected values equal to zero under Hq and s T C 1 s under 
Hi (see Fig. IA.ll) . For the same reason, the Pd is given by 



Pd=Q (V 1 (Pfa) - y/s T C- l 8^ . (A.13) 



Equation (IA.9I) can be written in the form 
T(x) = x T u > 7", 

with 

u = C^s. 



(A.14) 



(A.15) 



From this equation u can be thought as a linear filter of 
signal x. It is called matched filter (MF). 

A. 1.1. Some comments on the use of the matched filter in 
practical applications 

There are some important points to stress about the MF 
when used in practical applications, such as: 



T(x) is a sufficient statistic (jKav. I [l998). Loosely 
speaking, this means that T(x) is able to summarize 
all the relevant information in the data concerning the 
decision (IA.2I) . No other statistic can perform better; 



If the amplitude "a" of the source is unknown, 
Eq. (|A.9[) can be rewritten in the form 



then 



T{x) = x T C l g > 7" 



(A.16) 



with 7'" = i'/a = Q- 1 (P F A)v / 9 T C- 1 g. In other 
words, a statistic is obtained independent of "a" . 
As a consequence of the Neyman-Person theorem, in 
the case of unknown amplitude of the source, T(x) 
still maximizes Pd for a fixed Pfa- The only conse- 
quence is that Pd cannot be evaluated in advance. 
In principle this can be done a posteriori by using 
the maximum likelihood estimate of the amplitude, 
a = x T C~ 1 g/g T C~ 1 g. However, this is of little inter- 
est, since in real experiments one is typically interested 
in the detection of sources which have amplitudes 
characterized by a probability density function p(a). In 
this case, once Pfa is fixed to a value a and changing 
"a" across the domain of p(a ), the quantity 1 — Pd, 
with Pd given by Eq. (IA.13|) and s — ag, provides 
an estimate of the fraction of undetected sources as 
function of their amplitude; 

If s = Hs and x = Hx, then 



T{Hx) = xC s 

= x T H T H- T C- l H- l Hs = T{x), (A.17) 

with H any invertible linear operator (matrix). A use- 
ful consequence of this property is that if signal x is 
convolved with a function (e.g., the beam of an instru- 
ment), this operation does not modify the optimality 
of MF. In this case, the operator H transforms the co- 
variance matrix C into HCH T . This fact is useful in 
situations where more signals are available that are ob- 
tained with different point spread functions. 



A. 2. One-dimensional signal and multiple-frequency case 

In the context of CMB observations, the complexity in- 
creases since there are M signals x% = Sk + Jifc, Sfc = CLk9ki 
k = 1, 2, . . . , M, coming from the same sky area that are 
taken at different observing frequencies. Here, afc is the 
amplitude of the source at the A;th observing frequency, 
whereas g k is the corresponding spatial profile. For ease 
of notation, all the signals are assumed to have the same 
length N. In general, the amplitudes {a k } as well as the 
profiles {g k } are different for different k. However, if one 
sets 



n 



n ,n 



1 > ' l 2 > ' 



(A.18) 
(A.19) 
(A.20) 



it is possible to obtain a problem that is formally identi- 
cal to that treated in the previous section. Hence, the MF 
is still given by Eqs. (IA.14[) - (IA.15|) and is named multi- 
frequency matched filter (MMF). The only difference with 
the classic MF is that now C is a (NM) x (NM) block 
matrix with Toeplitz blocks (BTB): 



'Ml 



A I 



' MM 



(A.21) 



i.e., each of the Cy blocks is constituted by a JV x 
Toeplitz matrix. In particular, Cu = Fi[riinf] provides 
the autocovariance matrix of the ith noise, whereas C ij — 
EjVijiiJ], i 7^ j, the cross-covariance matrix between the 
ith and the jth ones. 

In spite of these similarities, when M > 1, there are 
additional difficulties: T(x) cannot be written in a form 
equivalent to Eq. (|A.16[) . In the case of unknown ampli- 
tudes {afc} T(x) cannot be computed. In other words, if 
the spectral characteristics of the radiation emitted by a 
source are not fixed, then the MMF is not applicable. 



A. 3. Extension to two-dimensional signals 

The extension of MF to the two-dimensional signals X and 
<S is conceptually trivial. In the case of one-dimensional 
case, setting 



s = VEC[5]; 
x = VEC[AT]; 
n = VEC \Af] , 



(A.22) 
(A.23) 
(A.24) 



formally the problem results in the same as given by 
Eq. (|A.2|) . In the case of multi- frequency observations, the 
situation is more complex since MF has to be applied to M 
signals at the same time. However, with the notation 

s = VEC [VEC[5i], VEC[5 2 ], . . . ,VEC[5 M ]] ; (A.25) 
x = VEC [VECLYi], VEC[X 2 ], . . .,VEC[X M }} ; (A.26) 
n = VEC [VEC[A/\], VEC[A/" 2 ], . . . , VEC[Af M ]} , (A.27) 

it is also possible to obtain a problem that is formally iden- 
tical to that given by Eq. Q 



The only difference is that 



6 VEC[JF] is the operator that transforms a matrix J- into a 
column array by stacking its columns one underneath the other. 
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now in Eq. (|A~21~j) . C is a (MN p ) x (MN P ) block matrix. If 
the signals {<Si} are two-dimensional N r x iV c maps, then 
each of the blocks is constituted by a (N r N c ) x (N r N c ) 
block Toeplitz with Toeplitz blocks (BTTB) matrix. In par- 
ticular, Cu provides the autocovariance matrix of the ith 
map, whereas Cy, i 7^ j, the cross-covariance matrix be- 
tween the ith and the jth maps. 

Especially for the multi-frequency case, the implementa- 
tion of the MF is not trivial. Even for moderate size signals, 
the matrix C becomes rapidly huge. As a consequence, it is 
necessary to implement numerical methods which are able 
to exploit the specific structure of C. Typically, they are 
based on Fourier approaches. This subject, however, is be- 
yond the aim of the present work. 
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Fig. A.l. Probability density function of the statistic T{x) under the hypothesis T-Lq (noise-only hypothesis) and Hi 
(signal-present hypothesis). The detection-threshold is given by 7". The probability of false alarm (Pfa) and the probability 
of detection (Pd) are shown in green and yellow colors, respectively. 



